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ABSTRACT 

We present three-dimensional, adaptive mesh simulations of dwarf galaxy out- 
flows driven by supersonic turbulence. Here we develop a subgrid model to track not 
only the thermal and bulk velocities of the gas, but also its turbulent velocities and 
length scales. This allows us to deposit energy from supernovae directly into super- 
sonic turbulence, which acts on scales much larger than a particle mean free path, 
but much smaller than resolved large-scale flows. Unlike previous approaches, we are 
able to simulate a starbursting galaxy modeled after NGC 1569, with realistic radia- 
tive cooling throughout the simulation. Pockets of hot, diffuse gas around individual 
OB associations sweep up thick shells of material that persist for long times due to 
the cooling instability. The overlapping of high-pressure, rarefied regions leads to a 
collective central outflow that escapes the galaxy by eating away at the exterior gas 
through turbulent mixing, rather than gathering it into a thin, unstable shell. Super- 
sonic, turbulent gas naturally avoids dense regions where turbulence decays quickly 
and cooling times are short, and this further enhances density contrasts throughout 
the galaxy- leading to a complex, chaotic distribution of bubbles, loops and filaments 
as observed in NGC 1569 and other outflowing starbursts. 
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1 INTRODUCTION 

Theoretical work has long shown that supernovae in low- 
mass galaxies should produce energetic outflows that heat 
and enrich their environments (Larson 1974; Dekel & Silk 
1986; Vader 1986). Since then such outflows have been ob- 
served in starbursting galaxies of all masses (e.g. Axon & 
Taylor 1978; Heckman 1990; Bomans et al. 1997; Martin 
1998; 1999; Heckman et al. 2000; Schwartz & Martin 2004; 
Veilleux et al. 2005) and at all cosmological epochs (e.g. 
Franx et al. 1997; Pettini et al. 1998; 2001; Frye, Broad- 
hurst, & Benitez 2002; Rupke et al. 2005). The existence 
of these ubiquitous galaxy outflows has several important 
implications for galaxy formation. The ejection of heavy el- 
ements has been invoked to explain the strong correlation 
between mass and metallicity observed in low-mass galax- 
ies (e.g. Dekel & Silk 1986; Richer & McCall 1995; Mateo 
1998; Thacker et al. 2002; Tremonti et al. 2004; Erb et al. 
2006; Kewley & Ellison 2008). The suppression of gas ac- 
cretion onto starbursting galaxies (Benson et al. 2003) and 
onto neighbouring density perturbations (Scannapieco et al. 
2000; Scannapieco et al. 2002) has been shown to be cru- 
cial to reconcile the small number of observed dwarf galax- 
ies with the large number of low mass dark-matter halos 



in the favored cosmological model (e.g. Somerville & Pri- 
mack 1999; Cole et al. 2000; Benson et al. 2003; Dekel & 
Woo 2003). The ratio of the baryonic mass to the gravi- 
tating mass of galaxies has been found to be several times 
less than the cosmic ratio (Hoekstra et al. 2005; Mandel- 
baum 2006), meaning that either the baryons never fell 
into galaxies or that powerful galactic winds removed them. 
And, widespread galaxy outflows have proven to be essen- 
tial to understanding the history of the intergalactic medium 
(IGM), which is observed to be widely enriched with heavy 
elements (Tytler et al. 1995; Songaila & Cowie 1996; Rauch 
et al. 1997; Chen et al. 2001; Simcoe et al. 2002; Schaye et al. 
2003; Aracil et al. 2004; Adelberger et al. 2005; Scannapieco 
et al. 2006), that drastically change its cooling properties 
(e.g. Sutherland & Dopita 1993; Wiersma et al. 2008). 

It has also become clear that most of the outflows that 
play a key role in each of these processes are driven by core- 
collapse SN ejecta and winds from massive stars. These cre- 
ate hot, metal-enriched bubbles that expand into the inter- 
stellar medium (ISM) and eventually break out of the host 
galaxy in the form of bipolar outflows that can reach ve- 
locities of hundreds of km/s (Veilleux et al. 2005). During 
this process the bubbles sweep up cooler ambient gas that 
can also be blown out of the galactic disk. In fact, all cur- 
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rent velocity measurements made in starburst galaxies are 
of the entrained cooler material, which is measured from 
UV or optical emission or via absorption lines (e.g. Rupke 
et al. 2005; Martin 2005). However, theoretical models pre- 
dict that around 90% of the energy and metal content of 
the winds exist in the hot (T ^ 10 6 K) phase, which has 
only recently been detected in X-ray emission (Strickland & 
Heckman 2007). 

The fact that the majority of the most important ma- 
terial remains invisible has led to a wide range of assump- 
tions regarding the efficiency with which starburst-driven 
winds can eject metal-enriched gas. Some authors claim that 
only winds from dwarf galaxies can reach the IGM (Ferrara 
& Tolstoy 2000) while other claim that winds can also es- 
cape from massive galaxies (Strickland et al. 2004). Com- 
pounding this observational uncertainty is the theoretical 
problem that the starbursting disk itself is not well mod- 
eled as a single-temperature medium in hydrostatic balance. 
Rather the gas is both constantly cooling and condensing 
into molecular clouds and being stirred and heated by ion- 
ization fronts, stellar winds, and supernovae (e.g. McKee & 
Ostriker 1977). This combination of extremely short cooling 
times and constant driving leads to a supersonic medium in 
which the turbulent velocities exceed the thermal velocities, 
and turbulent eddies act to support the disk even as they 
compress a fraction of gas, driving star formation (e.g. Mac 
Low & Klessen 2004). 

To approximate this configuration, previous simula- 
tions have been forced to both artificially increase the 
mass-averaged ISM temperature and suppresses its cooling 
through a variety of relatively crude techniques. The two- 
dimensional outflow simulations in Mac Low et al. (1989) 
and Mac Low & Ferrara (1999), for example, included an 
empirical heating function that was tuned to balance cool- 
ing in their initial configuration, but taken to be linearly 
proportional to the density, such that it was overwhelmed 
by cooling within dense regions that developed during the 
simulation. Strickland & Stevens (2000) following Tomisaka 
& Bregman (1993), Tenorio-Tagle & Munoz-Tunon (1998), 
and Suchkov et al. (1994) imposed a minimum temperature 
of 6.5 x 10 4 K in their two-dimensional simulations to ac- 
count for turbulent support of the disk. D'Ercole & Brighenti 
(1999) imposed a temperature floor in their simulations, 
equal to the 4.5 x 10 3 K temperature of the initial ISM, and 
Fujita et al. (2003; 2004) employed a similar hard cutoff at 
10 4 K. Mori et al. (2002) studied supernova feedback with 
full atomic cooling in a spherical "pregalactic" system with 
a 2 x 10 4 K virial temperature that was low-enough that 
catastrophic cooling of the inital was not a problem. Finally 
Cooper et al. (2008) carried out simulations with an inhomo- 
geneous ISM model made up of dense clouds of T ^ 3 x 10 4 K 
material surrounded by much more tenuous 5 x 10 6 if, "halo 
gas" such that very little gas in their simulations was located 
near the peak of the cooling function. 

While the extremely short cooling times within the ISM 
make it impossible to model supernova feedback by simply 
adding thermal energy to the gas, at the same time, the 
range of physical scales involved does not allow for the di- 
rect simulation of supernova remnants within a galaxy-scale 
simulation. Again, many approximations have been made to 
avoid this problem: ranging from temporarily lowering the 
densities of heated particles (Thacker & Couchman 2000), 



to delaying their cooling (Gerritsen & Icke 1997), to imple- 
menting momentum kicks rather than heating (Navarro & 
White 1993; Mihos & Hernquist 1994; Scannapieco et al. 
2001; Springel & Hernquist 2003). These theoretical uncer- 
tainties have also lead to suggestions that direct driving 
by supernovae and stellar winds (e.g. Silk 1997) may not 
be effective at all, but rather the primary driver might in- 
volve additional physics such as radiation pressure on dust 
(Thompson et al. 2005), or non-thermal pressure caused by 
cosmic rays (Socrates et al. 2008). In fact some of the scal- 
ing relations from such alternative models have been shown 
to reproduce many of the observed trends seen in the inter- 
galactic medium (Oppenheimer & Dave 2006). 

In this paper, we present simulations of starburst-driven 
outflows from a dwarf galaxy using an entirely new ap- 
proach. Building on a method developed by Dimonte & Tip- 
ton (2006), we develop a subgrid turbulence model that ac- 
counts both for the turbulent support of the disk and the 
extra turbulent energy input that drives a global outflow. 
Properly accounting for turbulence leads to a model for the 
galaxy that has a realistic temperature and pressure distri- 
bution, while at same time including cooling throughout the 
simulation. Furthermore, adding supernova input as turbu- 
lent energy lets us include this contribution without track- 
ing supernova remnants directly or adopting arbitrary ap- 
proximations to avoid the energy from being immediately 
radiated away. This then allows us to make observational 
predictions for the structure of the highly-disturbed ISM as 
well as assess the role of turbulent mixing on the escape 
fractions of gas, kinetic energy, and metals. 

The structure of the paper is as follows: In §2 we de- 
scribe our galaxy model, feedback model, and subgrid model 
for supersonic turbulence. Whenever possible we tune our 
parameters to approximate the starbursting dwarf galaxy 
NGC 1569, whose outflow has been well-studied at a va- 
riety of wavelengths. In §3 we present the results of our 
simulations, examine their dependencies on run parameters, 
and assess their observational consequences. Conclusions are 
given in §4. 



2 METHOD 

2.1 Simulation and Model Galaxy 

All simulations were performed with FLASH version 3.0, a 
multidimensional adaptive mesh refinement hydrodynamics 
code (Fryxell et al. 2000) that solves the Riemann prob- 
lem on a Cartesian grid using a directionally-split Piecewise- 
Parabolic Method (PPM) solver (Colella & Woodward 1984; 
Colella & Glaz 1985; Fryxell, Miiller, & Arnett 1989). In 
all runs we simulated a galaxy made up of a gas+stellar 
disk, contained within a dark matter halo. As a prototyp- 
ical dwarf starburst, we chose parameters to approximate 
the nearby galaxy NGC 1569, which has been well observed 
and analyzed in a wide variety of wavelengths (e.g. Reakes 
1980; Israel 1988; Gonzalez Delgado et al. 1997; Martin 1998; 
Greggio et al. 1998; Heckman et al. 2001; Stil & Israel 2002; 
Martin, Kobulnicky, & Heckman 2002) Consistent with ob- 
servations of this and other starbursting galaxies, we consid- 
ered a gas distribution with nearly exponential radial and 
vertical profiles (e.g. Barazza et al. 2006). However, as in 
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Roediger & Briiggen (2006) we softened the distribution in 
both directions, in order to prevent steep density gradients 
in the galactic plane and centre. This gave a gas distribution 
of 
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where (r, z) are the radius and distance from the plane 
in galactic cylindrical coordinates, M gas is the total gas 
mass, the vertical scale lengths are a gas and b gas , re- 
spectively, and £ = 0.9159 is Catalan's constant. For 
r It a gas and z > 6 gas , this density distribution 
converges towards the usual exponential disk p(r, z) — 
o„„-2 sa l — exp(-r/a gaa ) exp(-|z|/b gas ), and we fix M gas = 

gas "gas 

2 x 1O 8 M (Israel 1988), a gas = 0.7 kpc, and fo gas = 0.2 kpc 
(Reakes 1980). 

We did not calculate the self-gravity of the gas explic- 
itly, but rather as in Roediger & Briiggen (2006) we mod- 
eled the gravitational potential of the gas and stars as a 
Plummcr-Kuzmin disk (Miyamoto & Nagai 1975; Binney & 
Tremaine 1987), which approximates the exponential dis- 
tribution in disk galaxies, but is much more manageable 
analytically. The gravitational potential in this case is 
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where a disk = a gas , 6 disk = b gas , and M disk = 3 x 10 8 M Q 
are the radial scale length vertical scale length, and total 
mass, respectively, and G is the gravitational constant. We 
also added a second contribution to this potential due to the 
dark matter halo in which the galaxy is contained. In this 
case we assumed a Burkert (1995) model (see also Mori & 
Burkert 2000), which is given by 
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where R is the radius from the centre of the galaxy in 
spherical coordinates, Ro is the core radius of the halo, and 
p d ,o = 3.08 x 10~ 24 gcm~ 3 (_R /kpc)~ 2/3 is the central den- 
sity of the halo. For this potential, the maximum circular 
velocity of the halo as a function of Ro is 
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35 km/s as 



and for our model dwarf we assumed u c ,m ax 
in NGC 1569 (Martin 1998). 

Outside of the disk, the gas was assumed to consist of 
a uniform medium, with a mean density of p am bicnt = 10 -28 
g cm -3 , which is ~ 200 times the mean z — cosmological 
density. This gas was taken to be non-rotating, and in hy- 
drostatic balance with the assumed gravitational potential 
asymptotically approaching T a mbient ~ 2 x 10 4 -K" at large 
distances. Finally, we assumed that metals were smoothly 
distributed within the galaxy, such that all cells within the 
galaxy were initially at a fixed metallicity, Zinit = 0.25z7q 
(Gonzalez Delgado et al. 1997). For simplicity, we defined 
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Figure 1. Initial set-up for our galaxy outflow simulations. Top 
Left: Density in the midplane as a function of galactocentric ra- 
dius. Top Right: Density as a function of height at r = (solid 
line) and r = 2 kpc (dot-dashed line). Centre left: Total pres- 
sure as a function of radius in the midplane (solid line), and the 
component of pressure in the midplane arising from turbulence 
(dashed line). Centre right: Total pressure as a function of height 
at r = (solid line), and r = 2 (dot-dashed line) and the compo- 
nent of pressure at r = arising from turbulence (dashed line). 
Bottom left: Circular velocity as a function of radius in the mid- 
plane. Bottom right: Circular velocity as a function of height at 
r = 2. 



the boundary of the galaxy at p = 10p am bient, taking Z = 
0.01^0 outside of this region. 

All our simulations were performed in a three- 
dimensional 25 kpc x 25 kpc x 30 kpc region, with outflow 
boundaries on all sides. For our grid, we chose a block size 
of 8 3 zones and an unrefined root grid with 10 x 10 x 12 
blocks, for a native resolution of 313 pc. The refinement cri- 
teria were the standard density and pressure criteria, and 
for our fiducial runs we allowed for 3 levels of refinement be- 
yond the base grid, corresponding to an effective resolution 
of 640 x 640 x 728 zones, 39 parsecs on a side. 

Having fixed the density distribution in the disk, its 
pressure, turbulence, and temperature distribution were set 
such that: (i) hydrostatic balance was maintained in the di- 
rection perpendicular to the disk plane; and (ii) throughout 
the galaxy T ^ 10 4 K as atomic cooling quickly radiates 
thermal energy away above this temperature. In radial di- 
rection, the circular velocity was set so that the centrifugal 
force balanced the gravitational force and pressure gradi- 
ents. Figure [l] shows the resulting radial and vertical pro- 
files of density, pressure and circular velocity for our fiducial 
dwarf-starburst model. The choice of parameters for this 
model are summarized in Table [T] 
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Table 1. Fixed model parameters. 



Component 


Parameter 


Value 


gas 


Cigas 


0.7 kpc 




h 


0.2 kpc 




-^gas 


2 x 1O 8 M 




SFR 


0.25 M /yr 




^gas 


0.25 Zq 


gravitational 


a disk 


0.7 kpc 


potential 


£>disk 


0.2 kpc 




M diBk 


3 x 10 8 M Q 


DM halo 


rDM 


2 kpc 




Vc 


35 km/s 



2.2 Star Formation and Feedback 

For each galaxy, we computed a star-formation rate per unit 
area, Esfr from the surface density of gas, E ga s according 
to the Kennicutt-Schmidt law 
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which provides an accurate fit over » 6 orders of magnitude 
in S g as (Schmidt 1959; Kennicutt 1998). In this case the 
total star formation rate is given by 
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which we assumed to last for 50 Myrs in all of our simu- 
lations. For our choice of parameters this gave an overall 
star-formation rate of 0.25 Mq yr -1 , which is comparable 
to that observed in NGC 1569 (Israel 1998; Greggio et al. 
1998). This corresponds to a total mass of 1.3 x 1O 7 M stars 
formed over the course of each simulation. 

Each starburst was accompanied by energy and metal 
input from supernovae. One estimate of this contribution 
comes from comparing the cosmic SN rate per comoving 
Mpc 3 as measured by Dahlen et al. (2004) with the cos- 
mic star formation rate density as measured by Giavalisco 
et al. (2004), which gives a core-collapse supernova rate of 
(7.5±2.5) x 10~ 3 SNe per solar mass of stars formed (Scanna- 
pieco & Bildsten 2005). An alternative estimate is to count 
the number of M ^ 8M stars, by assuming a Salpeter ini- 
tial mass function with upper and lower mass cutoffs equal to 
120 Mq and O.1M , respectively (Scannapieco et al. 2002) 
which gives 1 SNe per 136 Mq of stars formed. We compro- 
mised between these two values and assumed that 1 SNe was 
generated per 150 Mq of stars formed, releasing an energy 
of 10 ergs. Furthermore, we assumed that a fraction /„ of 
this energy was instantaneously deposited into the galaxy, 
such that during the starburst the mechanical energy input 
was 



2.2 x 10 41 ergss _1 



SFR 



(7) 



,M Q yr- 

where in our fiducial models / w = 0.7. As described in more 
detail below, this energy was added as turbulent kinetic en- 
ergy, rather than as gas thermal energy, vastly changing the 
rate at which it was lost to cooling. 

As in a real galaxy, supernova energy was deposited 
into the gas stochastically, approximating the patchy distri- 
bution of OB associations within which stars and SNe are 



formed. In nearby galaxies, the luminosity function of OB 
associations is well-approximated by a power law of the form 



dNp-B 
dN 
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where Nob is the number of OB associations containing TV 
many OB stars, and /3 « 2 (McKee & Williams 1997; Oey & 
Clarke 1997). To approximate this distribution in our sim- 
ulations, we drew random numbers a* £ [0, 1] such that for 
each forming OB association i, the number of OB stars was 
given by 

-P) 



N t = + (1 - a,)^] 17 ' 1 " 



(9) 



where we took jV m i n = 20 and JV max = 1000. By drawing two 
other random numbers bi,d G [0, 1] we assigned each OB as- 
sociation a random azimuthal angle given by <f>i — 2na and 
a random radial position given by the transcendental equa- 
tion n = fjagas/1.5 where n = ln[(l + fi)/(l — &,)]. All 
associations were assumed to be centered around the mid- 
plane of the galaxy, and we paused for a time 150A r i/ SFR 
between bursts to maintain the overall star formation rate 
of the simulation. 

For each OB association, we injected f w W 51 Ni ergs of 
turbulent kinetic energy into the simulation in a region of 
size radius rbubble, which was at least the size of the re- 
gion containing twice the gas mass converted into stars, but 
no smaller than 60 pc so as to avoid extremely high pres- 
sure, largely-unresolved regions that excessively slow down 
the computational time step. The turbulent length scale for 
each OB association, discussed in more detail below, was 
taken to be rbubbic- Each SN was also taken to deposit ad- 
ditional gas and metals into this region. Here we assumed 
that the average ejected total mass of 8 Mq, per SN, 2 Mq 
of which is made up of heavy elements. This is consistent 
with the average stellar yields from a range of SNe II simu- 
lations (e.g. Maeder 1992; Wooseley & Weaver 1995; Arnett 
1996; Tsujimoto et al. 1995; Nagataki 1998), although there 
are significant theoretical uncertainties between various es- 
timates. 

Metallicity-dependent radiative cooling was calculated 
in the optically thin-limit, assuming local thermodynamic 
equilibrium: 



= -(1-50(1- 



0.25v t 2 /c 2 ) 



(10) 



where E coo \ is the radiated energy per unit mass, p is the 
density in the cell, m v is the proton mass, Y is the helium 
mass fraction, /j, the mean atomic mass, and A(T, Z) is the 
cooling rate as a function of temperature and metallicity, 
Vt is the turbulent velocity, discussed below, and c s is the 
local sound speed. Here we made use of the tables com- 
piled by Wiersma, Schaye, & Smith (2009) from the code 
CLOUDY (Ferland et al. 1998), making the simplifying ap- 
proximation that the abundance ratios of the metals both 
within the galaxy and ejected by the supernovae occured 
in the solar proportions. Furthermore, we also account for 
unresolved substrucutre in highly turbulent regions, assum- 
ing that within each cell the averaged density squared, (p 2 ) 
is given by the bulk density squared p 2 times an enchanc- 
ment factor that increases with the Mach number of super- 
sonic turbulence as 1 + 0.25(ut/c s ) 2 as measured by Padoan, 
Nordlund, & Jones (1997). 
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Because the disk in our simulations was initially sup- 
ported primarily by turbulent pressure rather than thermal 
pressure, and because the supernova energy was added to 
turbulence rather than to the thermal motions, no approxi- 
mate fixes to this equation were necessary. Instead, we were 
able to implement cooling in every cell in the simulation at 
every time step. 



2.3 Turbulence Modeling 

While the direct simulation of turbulence is extremely chal- 
lenging, computationally expensive, and dependent on res- 
olution (e.g. Glimm et al. 2001), its behavior can be ap- 
proximated to a good degree of accuracy by adopting a 
subgrid approach. Recently, Dimonte & Tipton (2006) de- 
scribed a subgrid model that is especially suited to cap- 
turing the buoyancy-driven turbulent evolution of subsonic 
bubbles (Scannapieco & Briiggen 2008; Briiggen & Scan- 
napieco 2009), which we have modified heavily to apply to 
supersonic galaxy-scale outflows. Other recent efforts at the 
subgrid modeling of turbulence on galactic and extragalac- 
tic scales are described in Maier et al. (2009) and Shen et 
al. (2009). 

Our model is based on the Navier-Stokes equations ex- 
tended to include a turbulent viscosity pr that depends on 
the eddy size L and kinetic energy per unit mass K. The 
interaction between the turbulence and the mean flow, is 
modeled by decomposing the flow into average components 
and fluctuating components, for example u to t = u + u' the 
sum of the mass-averaged mean-flow velocity and the fluc- 
tuating component of the velocity, and ptot = p + p' the sum 
of the mean-flow density and the fluctuating component of 
the density. The sum of these two components is substituted 
back into Navier-Stokes equations, which are then averaged 
to obtain separate evolutionary equations for the mean and 
fluctuating components. For compressible flows, the averages 
are weighted by the density such that 



ptotu' = 0, 



and 



Ptot Utot 



(11) 



(12) 



The resultant equations constitute an expansion about the 
mean flow that must be terminated with a simplifying set 
of closure assumptions. 

To leading order, the mean flow fluid equations in this 
case are given by 

dp dpUj 



dt dxj 

dpui dpinuj dP 
+ — - + 



= 0, 



dt 



dxi 



dpE dpEuj dPuj 
dt dxj dxj 



(13) 

0, (14) 

d fp^dPTX 
dxj \ N E dxj J 

+pEmcch + P-Ecool, (15) 



where t and x are time and position variables, p(x, t) is the 
average density field, «i(x, t) is the mass-averaged mean- 
flow velocity field in the i direction, P(x, i) is the total pres- 
sure, both turbulent and kinetic, and -E(x, t) is the mean 
internal energy per unit mass, also including both turbulent 



and thermal motions, and Ne = 1- Note that by modeling 
the impact of turbulence on the momentum equation by us- 
ing on ly a pressure, we are neglecting off-diagonal terms 
(dui/duj + duj/dui, where i ^ j.) This is because in the 
presence of shocks, such strain terms become unphysically 
large {e.g. Gauthier & Bonnet 1990; Klem 2004). Although 
these shortcomings may be overcome with appropriate use 
of limiters (e.g. Wilcox 1994; Sinha et al. 2003), we set this 
aside as a future refinement to the basic method presented 
here. 



In eq. (151, the first term on the right captures the ef- 



fects of turbulent mixing, which is modeled as a turbulent 
viscosity pt. the second, i? mC ch term is an explicit source 
term that is determined by the mechanical luminosity as 
given by eq. Q and the size of the region into which the en- 
ergy is being added, and the third, -B C ooi, term is given by eq. 
(10 1. Note that fit does not appear in the continuity equa- 



tion due to eq. (Ill and does not appear in the momentum 



equation due to eq. ( 12 \ 



In any type of subgrid model a number of fit parameters 
such as Ne arise, whose values are expected to be « 1, but 
must ultimately be fine tuned versus experiments to achieve 
the most accurate results. Here we take the same fit param- 
eters as used in Dimonte & Tipton (2006), and in Table 2 we 
summarize how they have been determined. However, there 
are several important differences in our approach in this pa- 
per, and thus it is likely our model will eventually be able to 
be further improved by re-adjusting these values to experi- 
ments and observations. Unlike in Scannapieco & Briiggen 
(2009), E is now the total internal energy, including both 
turbulent and thermal contributions. Likewise the pressure 
is the sum of both thermal and turbulent component, com- 
puted as P = „ 2p E, with m„ the mass of the proton and 
p the mean atomic weight of the gas. This redefinition al- 
lows us to apply the PPM solver to capture both the effects 
of turbulent and thermal pressures, yielding accurate solu- 
tions even in cases in which most of the internal energy is 
in the subgrid turbulent flow. An implicit and simplifying 
assumption associated with this approach is that turbulent 
and thermal velocities provide pressure support through the 
same equation of state. To track the metals ejected by su- 



pernovae, eqs. ( 13 1-( 15 1 are supplemented by a mass-fraction 
equation: 



dpF r t dpF r Uj 
~dt 



+ 



d.t 
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dxj 



Pt dF r 
N F dxj 



(16) 



where F r is the mass fraction of species r in a given zone, 
and Nf = 1.0 is a scale factor. 

The turbulence quantities that appear in these equa- 
tions are calculated from evolution equations for L, the scale 
of the largest turbulent eddies, and K, the turbulent kinetic 
energy. Simple equations for the evolution of these quantities 
are given by 



dpi 



+ 



dpLuj 
dxj 



d 
dxj 



pt dL 

NLdXj 



+ C c pL 



dui 
dxi 



and 
dpK 



+ 



dpKuj 
dxj 



o 

dx. 



( dK \ 



+ pE mcc h 



(17) 



(18) 



K_ dPuj 
' E dx. 



- pVtCn 



maxfVt-V^cO) 2 



where Vt = v2K is the average turbulent velocity, Vt,o is the 
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Table 2. Summary of coefficients in the turbulence model. Con- Table 3. Run Parameters 

stants that are fit to experiment appear with error bars, and in 

those cases we take the central value for this study. 

Run / w Resolution Subgrid Clumpy 

Parameter Value Effect Source Name (pc) 



N L 0.5 ±0.1 Diffusion of L Experimental Fit 7-4D 0.7 39 Yes No 

N E 1.0 ±0.2 Diffusion of E Self-similarity 4-4D 0.4 39 Yes No 

N F 1.0 ±0.2 Diffusion of Species Self-similarity 2.5-4D 0.25 39 Yes No 

N K 1.0 ±0.2 Diffusion of K Self-similarity 4-4N 0.4 39 No No 

C c 1/3 Compression of L Mass Conservation 7-3D 0.7 78 Yes No 

C D 1.25 ± 0.4 Drag term for K Experimental Fit 7-2D 0.7 156 Yes No 

7-4DC 0.7 39 Yes Yes 



average turbulent velocity at the start for the simulation, 
and Cc = 1/3 is given by mass conservation, Nl = 0.5, 
Nk = 1.0, and Cd = 1-25 are experimental fit constants. 
In the L equation the two terms on the right hand side rep- 
resent, respectively: turbulent diffusion and the growth of 
turbulent motions due to the expansion of the mean flow. In 
the K equation the three terms on the right hand side repre- 
sent, respectively: turbulent diffusion, the energy input from 
SNe, and a term that causes turbulence to decay away at a 
characteristic time scale of w L/(Vt — Vt,o) in the absence 
of external driving. This energy goes directly into gas heat- 
ing, because of energy conservation E — K ± ^thermal- The 
Vt — Vt,o in this term means that after the starburst episode, 
the gas remaining in the galaxy will eventually settle back 
into a turbulently-supported disk, rather than dissipate all 
its turbulence into heat. Finally, the turbulent viscosity was 
calculated as 



fit = pL max(Vj - V t ,o, 0), 



(19) 



where again the Vt — Vt,o term assures that in its initial 
configuration the galaxy remains static and does not diffuse 
into the intergalactic medium. Note that including the Vt,o 



term in eqs. ( 19 1 and ( 19 1 is akin to asumming a low level 



of persistent turbulence in addition to the much larger con- 
tribution by the SNe that are added to the galaxy explicitly 
during the simulation. 

Note that unlike in Dimonte & Tipton (2006) and our 
previous modeling, we do not include terms that attempt 
to track the growth of the Ray leigh- Taylor or Richtmyer- 
Meshkov instabilities on subgrid scales. This is because we 
are no longer working in the regime in which turbulence 
generated by these instabilities is dominant. Rather, we are 
interested in approximating the evolution of the ejecta asso- 
ciated with a random collection of SNe that move superson- 
ically within overpressured bubbles surrounding OB associ- 
ations. Our approach assumes that for each OB association, 
the initial maximum turbulent length scale is approximately 
equal to the radius of the overpressured region, and the ini- 
tial turbulent kinetic energy per unit mass is approximately 
equal to mechanical energy input from SNe divided by the 
total mass in the region. Although such overpressured re- 
gions will eventually expand and become Rayleigh- Taylor 
unstable, the length scale and velocities generated by this 
instability will initially be much smaller than the motions 
within the bubble. This means that the turbulent regions 
will not grow according to the linear equations of growth of 
the Rayleigh- Taylor or Richtmyer-Meshkov instabilities. In 
particular, one would not expect L and K in eqs. \17\ and 
(191, which represent SNe driven random motions, to grow 



as quickly as a t as they would in the Rayleigh- Taylor 
case. So our approach is to conservatively assume that L 
expands along with the mean flow, and that K is driven 
only by SNe input, and decays to thermal energy on the 
L/Vt timescale expected from a self-similar scaling analysis 
of incompressible turbulence (Kolmogorov 1941), which has 
also been confirmed in the compressible (magnetohydrody- 
namic) case (Stone et al. 1998; Mac Low et al. 1998; Padoan 
& Nordlund 1999). Thus only resolved Rayleigh- Taylor and 
Richtmyer-Meshkov instabilities are captured in our simula- 



tions, while eqs. (17 1 and (191 track the much more impor- 



tant subgrid turbulent motions driven by SNe. 

As in our previous modeling, our numerical implemen- 
tation of these equations was divided into three steps, which 
were carried out after the main hydro step in FLASH3, 
which advects all the variables above. In the first step, we 



implement the dui/dxi terms in eq. (171 explicitly. In the 
second step, we: (i) compute Vt as V2K, (ii) use a leapfrog 
technique to add the source term to Vt as Sk/pV, and then 
(iii) write V back to the K array as K = Vt/2. Finally, m 
the third step, we calculate the turbulent viscosity and use 



this to implement the diffusive mixing terms in eqs. ( 15 1-| 19 1 



explicitly. This final step requires us to impose an additional 
constraint on the minimum times step of dt ^ (A 2 pj pt)/i 
where (A) is the minimum of dx, dy, and dz in any given 
zone. This diffusive constraint must be satisfied for all zones 
in the simulation, but as discussed in SB08 this explicit ap- 
proach works well in concert with the AMR hydrodynamic 
solver, because as p t increases near the interface, density 
and pressure fluctuations are smoothed, allowing the code 
to derefine in these regions. Thus the diffusive time step re- 
mains greater than or comparable to the one required by the 
Courant condition for most of the evolution in our simula- 
tions. 

In all zones within the galaxy, K is initialized such that 
the galaxy is in hydrostatic balance, even though the initial 
thermal energy is chosen such that T ^ 10 4 everywhere. 
Outside the galaxy, the internal energy was taken to be 
solely thermal such that K was initialized to be to 10~ 10 iJ. 
Throughout the simulation, L was initialized to 2 parsec, 
l/100t/i of the initial gas scale height of the galaxy. 



3 RESULTS 

3.1 Fiducial Model 

In Figs. [2] - [4] we show results from our fiducial run, which 
has / w = 0.7 and 4 total refinement levels, including the base 
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Figure 2. Horizontal slices through the central 8x8 kpc in our fiducial simulation (7-4D) at four representative times: 10, 20, 30, & 
40 Myrs (arranged from left to right in each row). Top Row: Contours of log density from p = 10 -29 to 10 — 23 g cm -3 . Second Row: 
Contours of log temperature from T = 10 3 to 10 8 5 K. Third Row: Contours of log turbulent velocity, Vt, from 1 to 3000 km s _1 . Bottom 
Row: Contours of log turbulent length scale, L from 30 to 1000 parsccs. 



grid. A summary of all the runs carried out for this study is 
given m Table [3] Each run is labelled m-nD for cases with 
subgrid diffusive mixing and m-nN for cases without subgrid 
mixing, where m = 10 x / w and n is the number of refinement 
levels. Thus our fiducial run is referred to as 7-4D. 

Fig. [2] shows horizontal slices through the central 8x8 
kpc for this run at four representative times, and Fig. [3] 
shows vertical slices through the central 8 x 12 kpc. These 
plots contain numerous "superbubbles" driven by SNe from 
individual OB associations, as have been studied in several 
classic theoretical papers (Weaver et al. 1977; McCray & 
Snow 1979; Tomisaka & Ikeuchi 1986; Mac Low & McCray 
1988; Tenorio-Tagle & Bodenheimer 1988; Mac Low et al. 
1989). Note that the bubbles are somewhat larger at greater 
distances from the galactic centre as the ambient pressure is 
lower there. Regardless of galactocentric distance, however, 



the bubbles are significantly less empty than those described 
in the papers mentioned above. This primarily because we 
assume that SNe form and deposit energy into a region con- 
taining twice as much gas mass as is converted into stars. 
This means that for every 10 ergs added in SN driven tur- 
bulent energy, the affected region initially contains at least 
2x 15OM0 of gas. In this case, the increase in internal energy 
per unit mass is ~ / w 500 times that of the gas within the 
midplane, which corresponds to an initial turbulent velocity 
of « /w 1/2 500 km/s. 

Within these regions three time scales are almost equal: 
(i) the time scale for turbulent energy to dissipate into ther- 
mal energy, £diss ~ J*bubble/Vt; (u) the time scale for bub- 
ble expansion, t dyn ~ noubbie/c s , c ff ~ fbubble/Vt, where the 
effective sound speed within the heated regions is domi- 
nated by the turbulent motions, and (iii) the time scale 
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Figure 3. Vertical slices through the central 8 X 12 kpc in our fiducial simulation (7-4D) at four representative times of 10, 20, 30, & 
40 Myrs (arranged from left to right in each row). Top Row: Contours of log density from p = 10 -29 to 10 -23 g cm -3 . Second Row: 
Contours of log temperature from T = 10 3 to 10 8 ' 5 K. Third Row: Contours of log pressure from p = 10 — 17 ergs cm~ 3 to 10~ 10 ergs 
cm -3 . Bottom Row: Contours of log turbulent length scale, L from 30 to 1000 parsecs. 
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for mixing of the bubble interior with the surrounding mat- 
ter, U W rl ubhlc p/fl t ~ rb u bbl«/(VtrbubbJa) ~ nmbble/V*. 
Physically, this similarity of time scales occurs because the 
driving of gas into the exterior medium and the mixing of 
the exterior medium with the bubble interior both occur 
roughly at the time at which SNe kinetic energy is con- 
verted into thermal energy through a turblent cascade. Fur- 
thermore, all these time scales are short, on the order of 10 
pc / (/ w 1,/2 500kms~ 1 ) ~ 0.2 Myrs, which is much less than 
their interior cooling time and the dynamical time of the 
surrounding medium. 

Thus the bubbles quickly expand to the point that they 
are in pressure equilibrium with their surroundings, with 
only moderate heating of their interiors through turbulent 
decay, and moderate mixing of shells with the interior gas 
through turbulent "diffusion." As p oc n 5 / 3 oc £ -1 / 5 during 
this approximately adiabatic expansion, the result is the for- 



mation of w (/ w 500) 1/b x 60 pc « /w /6 200 pc superbubbles 
whose exterior shells are relatively thick, and whose interi- 
ors are roug hly (/ w 100) 3/5 « /w /5 40 times underdense and 
somewhat hotter and more turbulent than their surround- 
ings. As the turbulent length scale expands along with the 
flow, L also rises to ~ 200 pc within these regions, but re- 
mains well below the grid scale throughout the rest of the 
simulation. Many such regions can be seen in Figs. [2] and 
[3] particularly at later times and larger galactocentric radii. 
Note however, that turbulence will tend to be strongest in- 
side of the bubbles and outside of the swept up regions, as 
the time scale for turbulence to decay to thermal energy is 
much smaller at high densities. This can be seen from eqs. 



(171 and (19 1, which show that as pL and pK diffuse into 



denser regions the time scale for turbulent dissipation drops 
as tdocay oc L/V2K oc p~ 1 /y / p~ 1 oc p~ 1/2 . 

Furthermore the fact that atomic cooling is imple- 
mented throughout the simulation means the medium is un- 
stable to the cooling instability, and the dense shells around 
the bubbles will persist and grow over time. As described by 
Fall and Rees (1985) in the absence of thermal conduction 
and turbulence, dense clouds will become cooler and more 
compact over time if 



A(T P ,Z) A(Txsm,Z) 



± P 



T 2 
-'ism 



(20) 



where T p is the temperature of the perturbation, Tism is 
the temperature of the outside medium, and A(T, Z) is the 
cooling function. In the case of atomic cooling with Z ~ 
O.IZq, this is true whenever T p < Tism and T > 3 x 10 5 iC, 
which is clearly the case in the hot rarefied medium that 
develops within the galaxy. 

In our simulations turbulent pressure is also included, 



and eq. (201 is modified slightly to include the additional 



contribution to p dV work, yielding 

A(T P ,Z) A(T ISM ,Z) 
T${l + K p /E p y T? sm {l + K lsm /E lSM ) 2 ' K ' 

where K p /E p is the fraction the internal energy in turbu- 
lence within the perturbation, and Kism/Eism is the frac- 
tion of internal energy in turbulence in the exterior medium. 
As turbulence tends to persist longer in more rarefied me- 
dia, eq. (211 is even more easily satisfied than eq. (201. 



This means that unlike in simulations without cooling, the 
primary source of structure in our simulations is not the 



Rayleigh- Taylor instability that leads to the fragmentation 
of an initially smooth shell, but rather the stochastic nature 
of SN heating, enhanced by the cooling instability. Thus gas 
condenses into denser structures long after pressure equilib- 
rium has been achieved between the swept up gas and the 
exterior medium, as can be seen by comparing p and T to 
the pressure as in Figure [3] Note that a similar evolution 
of structure was also seen by Mori et al. (2002), who imple- 
mented cooling throughout their simulations of small "pre- 
galactic systems" and by Cooper et al. (2008), who imple- 
mented cooling in a larger galaxy, modeled as a hot rarefied 
medium surrounding distribution of cold, compact clouds. 

Fig. 2 also shows a large diffuse region that develops 
over time near the centre of the galaxy. Here the superbub- 
bles begin to overlap as star formation is strongly centrally 
concentrated due to both the overall radial profile of the 
gas and the E 1 ' 5 dependence of the Kennicutt-Schmidt law. 
Together, these overlaping outbursts open a rarefied region 
that expands gradually over time. As it grows, the collective 
outflow eats away at the exterior gas through turbulent mix- 
ing, rather than gathering it into a thin, fragile shell. As a 
result, the rarefied region drills its way almost directly verti- 
cally, following the path along which the minimum amount 
of material separates the bubble interior from the intergalac- 
tic medium. "Blow-out" occurs when the overpressured re- 
gion resulting from overlapping OB associations moves into 
the intergalactic medium, rather than when the material 
surrounding a single superbubble becomes Rayleigh- Taylor 
unstable. This occurs roughly at the time of bubble overlap, 
which can be estimated as 



overlap 



(R) m Abubbic(-R) 1 ^obE S f R (-R), 



(22) 



where Mob is the mass of a typical OB association and 
j4bubbic ~ 2MoBE~ a 1 s (/ w 500) 2 / 5 is the area of the disk cov- 
ered by a bubble after it expands to reach pressure equi- 
librium with its surroundings. From eqs. |5| and jTJ this 
gives 



tovcriap(-R) « 15 Myrs / w 2/5 exp ^- 



1.4kpc 



(23) 



This means that for any value of the wind efficiency, bubble 
overlap will occur much more quickly near the centre of the 
galaxy. The result is a strongly bipolar outflow of hot, diffuse 
gas with an opening angle that increases over time, which is 
often called the free wind. Within this region, gas densities 
are low, Vt and temperture are at their highest values, and 
the turbulent length scale, L, increases to values approach- 
ing a kpc as the gas rapidly expands above and below the 
disk. 

As discussed in Strickland et al. (2004), there is consid- 
erable observational evidence supporting the idea of super- 
winds being driven by the collective input of all the mas- 
sive stars near the central regions of the galaxy, rather than 
by the Rayleigh- Taylor break-up of individual superbubbles. 
Observed superwind pressure profiles, for example, demon- 
strate that mass and energy are ejected relatively smoothly 
over large regions (Heckman et al. 1990; Lehnert & Heckman 
1996), and the edges of well-resolved outflows match up well 
with the edges of starbursting regions (Strickland & Stevens 
2000; Strickland et al. 2000). We address the observational 
consequences of our results further in §3.3 below. 

Fig. [4] shows vertical contours of the late-time evolu- 
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Figure 4. Vortical slices through the central 20 X 30 kpc in our fiducial simulations at three representative times of 40, 50, & 60 Myrs 
(arranged from left to right in each row). Top Row: Contours of logp from p = 10 -29 to 10 -23 g cm -3 . Centre Row: Contours of total 
log energy density (kinetic+internal) from etot = 10 — 16 to 10 — 9 ergs cm -3 . Bottom Row: Contours of log total metal density from 

Pmctals = 10" 31 to 10~ 26 g cm" 3 . 



tion of the starburst. After blow-out, the collective outflow 
remains overpressured with respect to the surrounding inter- 
galactic medium, even at large distances. While the highly 
rarefied free wind does eventually collect up a denser shell of 
intergalactic gas during this expansion, mass loading arises 



mainly from dense clumps of interstellar material that are 
gradually mixed into the diffuse gas. Much of this material 
comes from the conical shear interface between the free wind 
and the surrounding galaxy, but there is also a contribution 
from clumps of gas being evaporated directly in the path 
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Figure 5. Rendered images of log density of metals from 10 -30 to 10 — 26 g cm -3 , separated into metals from the IGM (upper panels) 
and from the SNe that go off during the simulations (lower panels), at times of 20, 30, and 40 Myrs. The galaxy is being observed at 
an inclination of 30 degrees, and the galaxy disk visible in the upper panels has a radius of ss 6 kpc. Metals from SNe driving the 
outflow are ejected to large distances while the ISM is largely retained by the galaxy. The plotting artifacts at large radii in the top 
panels illustrate the 8 3 cell "blocks" in our adaptive mesh approach, which has a resolution of 39 parsecs in the central regions where 
superbubbles develop, but only 312 parsec at large radii, where density and pressure contrasts are small. 



of the outflow. Together these mixed, entrained components 
provide most of the observational constraints on starburst- 
driven winds. 

Metals can escape from the galaxy either by being en- 
trained in the wind or by being directly ejected in the SNe 
remnants driving the outflow. To distinguish between these 
two contributions, in Fig, [5] we show rendered images of log 
Pmetais separated into the component arising from the metals 
initially in the ISM, and metals that arise from SNe occur- 
ring during the simulation. Here we see that despite the fact 
that turbulence is driving the outflows in our galaxies, the 
mixing of supernova ejecta into the ISM is minimal, and each 
of the two metal components evolves completely differently. 
Consistent with the blow-out picture of outflow generation, 
the ISM metals are swept up into dense shells of gas that re- 
main bound to the galaxy. On the other hand, SNe ejecta are 
only weakly mixed into the shells, such that the metallicity 
changes very little within dense regions during the starburst, 
as observed in the HII regions in NGC 1569 (Kobulnicky & 



Skillman 1997). Instead SNe metals are mostly found within 
the rarefied high-pressure regions, and are able to escape to 
large distance following bubble overlap. 

To quantify the kinematics of the ejected gas further, 
we plot in Fig. [6] the evolution of the ejected mass, energy, 
momentum, and metals. In the upper panel of this figure, 
we show the gas mass ejected by the galaxy as a function 
of time, compared to the total initial gas mass. Here we de- 
fine escaping gas as material that is at least a scale height 
from the midplane of the galaxy and traveling outwards such 
that the component of velocity in the direction of the grav- 
itational acceleration exceeds the local escape velocity. 

The wind is able to effectively blow out ~ 3 x lO 6 A/0 
of gas from the dwarf starburst, which is comparable to the 
total mass of stars formed, as observed in local starbursts 
(Martin 1999). However, the ejected mass is relatively small 
as compared to the total 2 x 10 s Mq gas mass of the galaxy, 
and thus the majority of the ISM remains bound. This is 
true even though the total kinetic energy from SNe added to 
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Figure 6. Evolution of ejected fractions in our fiducial simulation 
(7-4D). Top Panel: Evolution of the ejected gas mass (solid line) 
as compared to the total gas mass of 2 X 10 8 Mq (thick horizontal 
line). Second Panel: Evolution of the ejected kinetic energy (solid 
line) and ejected internal energy (dashed line) as compared to 
the energy input from SNe (dotted line) and M ga3 v^ BC /2 (thick 
horizontal line). The small decrease in ejected KE at 55-60 Myrs 
is due to the some of the fastest moving material moving out of 
the simulation volume. Third Panel: Evolution of the total ejected 
momentum (solid line) as compared to Afg as t) csc (thick horizontal 
line). Bottom Panel: Evolution of the ejected mass of metals (solid 
line) and metals coming purely from SNe going off during the 
simulation (dashed line), versus the total mass of metals added 
to the simulation (dotted lines). 



the simulation, shown in the second panel of Fig. [6] exceeds 
the binding energy of the galaxy. Rather the majority of 
the energy deposited near the galaxy centre is carried away 
vertically by the outflowing wind, while most of the energy 
in superbubbles at larger radii decays to thermal energy and 
is radiated away. In fact, as discussed in Mac Low & Ferrara 
(1999), "blow-away" of the ISM only occurs when the lateral 
walls of the central outflow are accelerated so quickly that at 
the end of the blow-out phase they are moving with enough 
momentum to sweep out the rest of the ISM gas radially. 
This requires significantly more energy than deposited here, 
although the condition for blow-away is far less restrictive 
for relatively round and puffed up galaxies such as ours than 



it is for larger and less turbulent disks (Mac Low & Ferrara 
1999). 

In the third panel of Fig. [6] we compare the total mo- 
mentum of the ejected gas, which is at least an order of 
magnitude smaller than the total galaxy gas mass times the 
escape velocity. At late times K -Ejected /-Pejected ~ 200 km/s 
indicating that the majority of the outflow escapes at a few 
times the escape velocity, although again most of this mass 
is in the form of the largely invisible free wind rather than 
the much better observed entrained material. 

In the bottom panel of this figure we compare the metal 
mass produced during the starburst to the total ejected 
metal mass and the ejected metal mass coming from SNe 
that occur during the simulation. As seen in Fig.[5]the ma- 
jority of metals ejected from the galaxy come from the super- 
novae that are driving the wind itself, while only a small frac- 
tion comes from metals already contained within the ISM. 
However, most of the metals produced by the starburst are 
retained by the galaxy, confined to hot regions that mix into 
the dense ISM on time scales much longer than the starburst 
itself. 



3.2 Parameter Dependencies 

3.2.1 Wind Efficiency 

In Figs. [7] and [§] we study the impact of reducing the effi- 
ciency with which SNe energy is converted into turbulent 
energy (/ w ) and the efficiency of turbulent mixing (C M ). In 
the top two rows of Fig. [7] we show vertical slices of density 
for runs 4-4D and 2.5-4D, for which / w has been reduced 
to 0.4 and 0.25 respectively. Here we see that varying / w 
leads to drastic differences in outflow strengths and mor- 
phologies. The galaxy in our fiducial (7-4D) run achieves 
blow-out within 30 Myrs from the start of the simulation, 
rapidly expanding to fill most of the simulation volume by 50 
Myrs, but reducing / w delays blow-out to « 40-50 Myrs in 
the / w = 0.4 and / w = 0.25 cases, and substantially reduces 
the volume enriched by the outflowing gas. This can be un- 
derstood in the context of a collective wind that arises from 



overlapping superbubbles. From eq. (231, we see directly 
that at a fixed radius, overlap occurs much later in models 
with low / w as the individual superbubbles are much more 
compact. Note also that the opening angle of the outflow is 
reduced as / w decreases, and the central collective outflow 
is only able to punch its way through a smaller region of the 



ISM. Again this is consistent with eq. ( 23 ) which shows that 
for a given time, it is only the most central regions that can 
achieve superbubble overlap, followed by "blow-out" of the 
gas almost vertically. Thus in the / w = 0.25 case, not only 
does a small fraction of the hot gas escape from the galaxy, 
this is concentrated into a plume traveling almost directly 
perpendicular to the disk. 

In the lower panel of this figure, we examine the effect of 
reducing turbulent mixing by setting / w = 0.4 and C M = 
in a case we label as 4-4N, where N indicates that no subgrid 
diffusive mixing has been modeled. In this run, the mixing 
between bubble interiors and the shells is severely reduced, 
which in turn reduces radiative losses, and results in some- 
what larger superbubbles. In general, this run behaves sim- 
ilarly to our fiducial case in which we assume more efficient 
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Figure 7. Vertical slices of log p from p = 10 29 to 10 23 g cm 3 , through the central 20 X 30 kpc in our comparison simulations (4-4D, 
top, 2.5-4D middle, and 4-4N bottom). From left to right t = 30, 40, 50, & 60 Myrs, respectively. 



SNe driving of turbulence, but mix a significant fraction of 
the energy into dense, rapidly radiating regions. 

In Fig. [8] we plot the ejected gas mass, energy, momen- 
tum, and metal mass for each of the runs shown in Fig. [7] As 
apparent from the density slices, the ejected mass and mo- 
mentum depend extremely sensitively on / w . After 60 Myr 
the galaxy with / w = 0.4 has ejected w 10% as much mass 
as in the fiducial / w = 0.7 case, while the / w = 0.25 galaxy 
has ejected only « 3%. The differences in ejected energy 
and momentum between the runs are even more dramatic. 
The energy ejected in the / w = 0.25 run, for example, is less 
than 1% that of the fiducial, / w = 0.7 run. This large dif- 
ference also translates into a large difference in the ejected 
metal fraction, which is almost negligible in the / w = 0.25 



case. Finally, in the no mixing case with / w = 0.4, ejection 
of mass, energy, momentum, and metals are all significantly 
enhanced, reaching values similar those in our fiducial run. 



3.2.2 Resolution 

Next we carried out a test of the impact of resolution on 
our results. Leaving the base grid spacing fixed at 313 pc, 
we resimulated our fiducial / w = 0.7 galaxy with different 
maximum levels of refinement. In the first of these runs, 
labeled 7-3D, we only allowed 2 levels of refinement above 
the base (level 1) grid, for an effective resolution of 78 pc. In 
the second run, labeled 7-2D, we allowed only 1 additional 
level, for an effective resolution of 156 pc. In this very low- 
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Figure 8. Evolution of the gas mass (top), energy (2nd row), momentum (third row) and metal mass (bottom) for runs 4-4D (left 
column), 2.5-4D (right column), and 4-4N (left column). Lines are as in Fig. [6] 



resolution case, the minimum value of rbubbie was taken to 
be 120 pc instead of 60 pc, so that each bubble region would 
encompass more than a single zone. 

Vertical slices at various times through the simulation 
volumes from these runs are shown in Fig. [9] Note when 
comparing these runs that the turbulent length scale can 
grow larger than the grid scale, as shown in Figures [2] and [3j 
and thus turbulent diffusion can smooth features on scales 
larger than the effective resolution of each of the runs. With 
this limitation in mind, we see that in general, both low- 
resolution runs display the same evolution as in the fiducial 
run. In all cases, a low density cavity builds up near the cen- 
tre of the galaxy, pushing its way through the lowest-density 
regions until it makes its way into the surrounding inter- 
galactic medium. While blow-out occurs at slightly different 
times in each of the runs, it is always an abrupt transition 
that is quickly followed by gas ejection out to very large 
distances. At late times, in all runs, the base of the outflow 
widens and the flow becomes less collimated as bubbles over- 
lap at larger radii and significant mass entrainment occurs 
at the interface between the free wind and the surrounding 
galaxy. However, it is only in the highest resolution run that 
clouds of dense ISM are resolved within the outflow itself. 

In Fig.[l0]we compare the evolution of the ejected quan- 
tities as a function of resolution. At early times, before the 
initial blow-out occurs, there are notable difference between 
the runs, and in general, the higher resolution cases achieve 
higher ejected fractions earlier. As the simulations progress, 



however, the evolution becomes much more similar between 
the runs, and by the final time of 60 Myrs, all ejected quan- 
tities are consistent to within less than a factor of two. 
This is true even though the maximum volume resolution, 
and hence the mass resolution, varies by a factor of 2 6 be- 
tween these runs. This means that our implementation of 
supersonic turbulence-driven outflows is only weakly depen- 
dent on the extent to which turbulence is directly resolved, 
as opposed to approximated by subgrid modeling. Instead, 
the qualitative and quantitive evolution of the starburst is 
largely independent of resolution. 

3.2.3 ISM Structure 

Next, we examined the effects of pre-existing ISM structure 
on the evolution of the outflow. As it lies outside of the 
scope of this paper, we did not attempt to model a realistic 
density distribution of molecular clouds. Rather we simply 
altered the gas distribution so as to add a regular series 
of dense « 0.5 kpc knots, which interact with the outflow 
as it develops. Specifically, we adopted a perturbed density 
distribution within the galaxy, given by: 



Pperturbed(£, y, z) = 



pavcragc 

(x, y, z)[l - 0.8 x 
cos(7ra;/A) 2 cos{-k V /\) 2 ] 2cos{7Tz/x \ (24) 



where A = 0.5 kpc and ^average is the smooth distribution 
given by eq. |T]). At the same time we rescaled the temper- 
ature and turbulent kinetic energy per unit mass through- 
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Figure 9. Vertical slices through the central 8 X 12 kpc in simulations with varying resolutions at four representative times of 10, 20, 
30, & 40 Myrs (arranged from left to right in each row). The lowest resolution run, 7-2D, is shown in the top row, the intermediate 
resolution run, 7-3D, in shown in the centre row, and the highest resolution, fiducial run, 7-4D is shown in the bottom row. All panels 
give contours of logp from p = 10 -29 to 10 — 23 g cm -3 . 



out the galaxy as [ppcrturbcd/pavcragc] -1 such that the over- 
all pressure profile remained the same as in the fiducial run, 
with each of the dense clouds in pressure equilibrium with its 
surroundings. All other parameters for this run were iden- 
tical to the fiducial 7-4D run, and we refer to it as 7-DC, 
where the C indicates the presence of a clumpy ISM. 

In Fig. [ITJ we show vertical slices of the density, tem- 
perature, and pressure from this simulation at several repre- 
sentative times. These illustrate the strong tendency for the 
outflow to avoid dense pockets of gas. Again, this is both 
because at hig densities radiative cooling, which goes as p 2 , 



is much more efficient and because the time scale for turbu- 
lence to decay to thermal energy is much shorter. 

As in the fiducial run, superbubbles in the 7-4DC run 
begin to overlap into a large rarefied region near the centre of 
the galaxy, but the asymmetry in this run is much stronger 
than in the fiducial case. Here, the gas moves primarily later- 
ally in the plane of the disk, so as to avoid passing directly 
through the pair dense clumps located directly above and 
below the centre. This is an extreme example of the cooling 
instability, which preserves dense clumps even in the direct 
path of the hot, high-pressure gas. The hot region pushes 
its way around these clumps, forming two distinct chim- 
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Figure 10. Comparison of the evolution obtained at varying maximum levels of refinement . As in Fig. [6] from top to bottom the panels 
show the evolution of the gas mass (top), energy (2nd row), momentum (third row), and metal mass (bottom). From left to right the 
models correspond to 1 levels of refinement beyond the base grid (7-2D) yielding 156 pc resolution, 2 levels of refinement beyond the 
base grid (7-3D; 78 pc resolution), and the fiducial case (7-4D; 39 pc resolution). Lines are as in Figs.|6]&[8] 



neys of hot material that punch out into the intergalactic 
medium on either side of the central axis. The free wind 
then streams outwards on both sides of dense regions, en- 
veloping this cold gas and eventually entraining it into the 
outflow. Only at late times does the gas within the clumps 
eventually begin to move out of the galaxy, as it gradually 
shears and mixes into the wind much like the ISM along the 
edges of the outflow. 

Note that this evolution is markedly different than what 
would occur if SN heating were modeled purely as thermal 
energy input and radiative cooling were neglected. In this 
case, the shocks from the developing outflow would raise 
the pressure in the cold clouds, causing them to expand and 
smoothing out the overall density distribution. A tendency 
to preserve a clumpy medium near the base of the wind, as 
observed in NGC 1569 (Hunter, Hawley & Gallagher 1993; 
Tomita, Ohta & Saito 1994; Heckman et al. 1995; Martin 
1998; Westmoquette, Smith, & Gallagher 2008), is one of 
the hallmarks of outflows driven by supersonic turbulence. 



3.3 Observational Consequences 

Finally, we carried out a preliminary comparison of our mod- 
els with observations of NGC 1569. Here we focus on Ha 
imaging as an optical tracer of warm, dense, ionized gas and 
X-ray imaging as a tracer of the hot, rarefied gas. 



3.3.1 Ha Emission 

As a tracer of the shocked and ionized interstellar medium 
we computed the Ha emissivity of the gas, )Ha, in each 
grid zone throughout the simulation. Note that for reasons 
of scope we neglect photoionization for this calculation, al- 
though this can be very important source in the presence 
of large numbers of O and B stars as would naturally oc- 
cur during a starbursts. For this reason the images below 
should be considered as only tracing the rough morphol- 
ogy of the ionized gas, while more detailed studies, such as 
comparisons between the velocity structure of Ha emitting 
gas and the coldest outflowing gas as measured though Nal 
absorption {e.g. Heckman et al. 2000; Fujita et al. 2009), 
will require more complete calculations including radiative 
transfer effects. 

With these limitations in mind Ha can be calculated as 
function of temperature and density from the H/3 emissivity 
of the gas from the theoretical fit (Ferland 1980) 



n c n p 



2.53 x 10" 
1.12 x 10" 



-o.833 erggg -i f or j; <; 2.6 x 10 4 iv" 
T^ergss -1 for T e > 2.6 x 10 4 A', 

(25) 

combined with an assumed Balmer decrement, jir a /jff/3- 
For simplicity, we fixed ]Ha/]Hp = 2.9, ignoring the small 
« 10% changes that occur over the range of temperatures 
encountered in starbursts. Again neglecting any photoioniz- 
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Figure 11. Vertical slices through the central 8 X 12 kpc in our perturbed simulations (7-4DC) at four representative times of 10, 20, 
30, & 40 Myrs (arranged from left to right in each row). Top Row: Contours of logp from p = 10 — 29 to 10 — 23 g cm -3 . Centre Row: 
Contours of logT from T = 10 3 to 10 8,5 K. Bottom Row: Contours of logp from p = 10~ 17 ergs cm -3 to 10 — 10 ergs cm -3 . 



ing background and assuming local thermodynamic equilib- 
rium, we calculated n e = n p directly from the Saha equation 
as 

ex P I — — I i ( 26 ) 



n h 



( 2Tim e kT 

{—IT 



kT 



where nu is the neutral hydrogen density, m e is the mass 
of the electron, and h is Planck's constant. To account for 
unresolved gas inhomogeneities we assummed this emission 
was enhanced by a factor of 1 + 0.25(Vt/c s ) 2 , as in eq. (10 1. 
Finally we projected the total Ha emissivity in the z and 
and x directions, to produce vertical and horizontal surface 
brightness maps of our fiducial (7-4D) starbursting galaxy. 



Plots of the logarithm of the surface brightness at times 
of 20, 30, and 40 Myrs are given in Fig. [12] As in NGC 1569, 
the simulated vertical Ha profiles shown in the upper panels 
of this figure exhibit a complex and chaotic structure, which 
is is brightest near the plane of the galaxy. Also as seen in 
NGC 1569, which is observed almost edge-on, our images 
display bubbles and loops of strong Ha emission, which cor- 
respond to the shells of material swept up by the superbub- 
bles generated by each OB association {e.g. Hunter, Hawley, 
& Gallagher 1993; Martin 1998; Westmoquette, Smith, & 
Gallagher 2008). 

Furthermore, and especially at late times, outgoing fil- 
aments of heated gas are apparent near the central axis 



© 2009 RAS, MNRAS 000, [T]-?? 



18 E. Scannapieco & M. Briiggen 





-1 o 1 
x/kpc 



Figure 12. H Q luminosity (arbitrary scale) of our fiducial simulations (7-4D) at three representative times of 20, 30, & 40 Myrs (arranged 
from left to right in each row). Top Row: Central 6x6 kpc 2 , projected in the z direction. Bottom Row: Central 6x6 kpc 2 , projected in 
the y direction. 



(Tomita, Ohta & Saito 1994; Heckman et al. 1995; Mar- 
tin 1998). This is also consistent with observations, and a 
comparison of these plots with the density contours in Fig. 
[3] shows that these features arise largely from ISM material 
that is being entrained by the hot wind, either because the 
dense gas is on the edges of the blow-out region, or because 
it is contained in cold clumps that are left directly within the 
path of the collective central outflow. Again, these clumps 
form naturally in our models from the material swept up be- 
tween superbubbles, and they are enhanced by the cooling 
instability and the tendency for supersonic turbulent gas to 
avoid density perturbations rather than disrupt them. 



In the lower panels of Fig. [12] we show horizontal, z- 
projected views of the Ha distribution in our simulations. 
While these cannot be directly compared to observations of 
NGC 1569, they nevertheless serve to illustrate the cellu- 
lar nature of the warm, ionized gas, which becomes more 
and more concentrated into the shells of superbubbles as 
time progresses. Note also the dark gap in the centre of the 
images, which grows over time. This is the signature of blow- 
out, which allows us to peer straight through the centre of 
the galaxy from this vantage point, unimpeded by the warm 
dense medium that previously kept the X-ray emitting cen- 
tral gas confined to the central starburst. 



3.3.2 X-ray Emission 

Next we turn our attention to the X-ray properties of our 
simulated galaxy. Such observations provide the most di- 
rect picture of outflowing starbursts as they reveal diffuse 
emission not only from the disk and the halo, but from the 
hottest and highest pressure regions of the flow. Martin et 
al. (2002), for example, have taken Chandra observations of 
NGC 1569 and found that its X-ray luminosity is dominated 
by diffuse, thermal emission from the disk (0.7 keV) and a 
bipolar 0.3 keV halo. After subtracting hard point sources, 
they found a luminosity from the thermal component in the 
band of 0.3-6 keV of « 8 x 10 38 erg s -1 . Moreover, Chandra 
observations of galactic winds have shown that the X-rays 
are often spatially-correlated with Ha emission (Cecil et al. 
2002; Strickland et al. 2004; Grimes et al. 2005). A promi- 
nent example is the edge-on spiral NGC 3079 (Cecil et al. 
2002), which shows towers of intertwined Ha filaments that 
line its central outflow and emit in the X-ray, as discussed 
by Strickland et al. (2002). 

From a theory point of view, Suchkov et al. (1994) com- 
puted the X-ray emission from two-dimensional simulations 
of galactic superwinds that interacted with a two-component 
ISM. They found that the bulk of the soft X-ray emission 
originated from shocked material in the disk and halo, and 
the bulk of the hard X-ray emission arose from the free wind 
itself. Consequently, they concluded that soft X-ray spectra 
need not show abundances enhanced in metals. Strickland 
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Figure 13. False colour maps of the logarithmic surface brightness in the soft (0.5-2 keV; red) and hard (2-10 keV; blue) X-ray bands 
in our fiducial simulation (7-4D) at 30 (left panels) & 40 Myrs (right panels). The top row shows a projection in the y-direction and the 
bottom row shows a projection in the z-direction, spanning the central 12 X 12 kpc. 



& Stevens (2000) also performed two-dimensional starburst 
simulations, computed the X-ray emission, and found that 
most of the soft X-rays come from shock-heated ambient gas 
and from the interface between the hot gas and the ISM. 
However, the symmetry of these simulations did not allow 
them to form filamentary structures to the same degree as in 
three-dimensional models, and this limited their predictions. 

For this study, we computed the soft (0.5-2 keV) and 
hard (2-10 keV) X-ray emissivities of our fiducial simula- 
tion using the ATOMDB cod^J which includes the Astro- 
physical Plasma Emission Database (APED) and the spec- 
tral models output from the Astrophysical Plasma Emission 
Code (APEC). The APED files contain information such 



http:/ /cxc. harvard.edu/atomdb/ 



as wavelengths, radiative transition rates, and electron col- 
lisional excitation rate coefficients, and APEC uses these 
data to calculate model spectra, for optically-thin plasmas 
in collisional ionization equilibrium. We also included sub- 



grid enhancement as in eq. ( 10 1 and projected these metal- 



dependent emissivities in each grid cell to give the surface 
brightness maps shown in Fig. |13| which shows the distribu- 
tion at 30 and 40 Myrs. In these composite colour images, 
the red colours correspond to the soft band, the blue colours 
correspond to the hard band, and both colour scales cover 
six orders of magnitude. 

This figure shows that most of the soft X-ray emis- 
sion comes from the the winds that are blown out of the 
disk. Also, we find towers of X-ray emission that rise more 
than 4 kpc above the disk and follow the channels through 
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which the hot wind escapes (cf. Fig.[3j). The hard X -rays, on 
the other hand, come mainly from the starbursting bubbles, 
most of which are in the centre of the galaxy. Comparing 
Figs. |12| and |13| it appears that the sites of X-ray emission 
are correlated, similar to the observations of NGC 3079. This 
is because the outflow is directly responsible for accelerating 
filaments of cold gas out of the disk. 

The X-ray emission is very temperature dependent and 
thus the emission changes over time. At 30 Myrs the soft 
X-ray emission dominates the emission from the wind and 
the regions that are emitting hard X-rays are confined to 
the starbursting bubbles. At this time, the overall X-ray lu- 
minosity is ~ 6 x 10 38 erg s _1 in the 0.5-2 keV band and 
« 3 x 10 35 erg s" 1 in the 2-10 keV band. The mean X- 
ray weighted temperature of the soft X-ray emitting gas is 
1.6 x 10 e K and of the hard band 7.6 x 10 7 K. The soft 
X-ray emission is so luminous because most of the emission 
measure occurs at temperatures at which emission is strong. 
Thus the X-ray luminosity is very close to the value mea- 
sured for NGC 1569. 

At 40 Myrs, the wind has moved much of the gas to 
higher temperatures with the effect that the emissivity in 
the soft band drops, while the emissivity in the hard band 
remains roughly constant. The overall X-ray luminosity is 
k 9 x 10 36 erg s _1 in the 0.5-2 keV band and « 4.5 x 10 35 
erg s _1 in the 2-10 keV band. The mean X-ray weighted 
temperature of the soft X-ray emitting gas is 5.6 x 10 6 K 
and of the hard band 3.7 x 10 7 K. At 40 Myrs, emission from 
the free wind can been seen in the hard X-ray band, moving 
out to large distances and escaping from the potential well 
of the galaxy (seen as blue patches in Fig. |13[ ) . 

However, one should note that we do not include a hot 
X-ray emitting halo, and have calculated only the thermal 
emission from the shocked disk and the supernovae. Espe- 
cially in the hard X-ray band non-thermal emission becomes 
important. Real spectra show a combination of thermal 
emission and non-thermal components due to X-ray bina- 
ries, young supernova remnants, low-luminosity AGN, and 
Compton scattering of relativistic electrons by the ambient 
far-IR and CMB radiation fields (see e.g. Persic & Rephaeli 
2002). Therefore, a detailed comparison to observations is 
not straightfoward, especially in the hard band, and for rea- 
sons of scope we postpone this to a future publication. 

Finally, we note that a major driver of future X-ray 
instruments such as the X-ray Microcalorimeter System on 
board of the International X-ray Observatory (IXO^ is to 
take high-resolution spectra of galactic winds in order to 
measure their velocity, temperature and abundance struc- 
ture. While models such as ours that are able to capture 
turbulent broadening of lines will be essential for this mis- 
sion, we do not attempt to show spectra in this study as 
a thorough simulation of the X-ray emission would need to 
include a number of additional important processes such as 
absorption. Again, this is left for future work. 



2 http://ixo.gsfc.nasa.gov/ 



4 SUMMARY AND CONCLUSIONS 

Starburst-driven outflows play a key role in structure forma- 
tion, impacting issues ranging from the gas and metallicity 
evolution of galaxies to the chemical and thermal history 
of the intergalactic medium. While much observational and 
theoretical progress had been made in understanding these 
objects, formidable challenges remain. Many local outflow- 
ing starbursts have been observed in great detail, but the 
phase that contains 90% of the ejected energy and metals has 
gone largely undetected. Many theoretical studies of galaxy 
outflows have been conducted, but these have been faced 
with large uncertainties arising from rapid radiative cool- 
ing and a complex turbulent gas distribution that contains 
structures over a wide range of scales. In fact, the medium 
within starbursting galaxies is disturbed so strongly, and the 
cooling times within the gas are so short, that the turbulent 
velocities far exceed the thermal velocities. 

Here we have explored a completely new theoretical ap- 
proach that addresses these issues by tracking not only the 
thermal and bulk velocities of the gas, but also its turbu- 
lent velocities, pressures, and length scales. By adding an 
intermediate class of gas motions that operates on scales 
much larger than the particle mean free path, but much 
smaller than the resolved motions in the simulation, we are 
able to carry out starburst simulations that overcome many 
of the problems seen in previous models. In particular, our 
approach allows us for the first time to model starbursting 
galaxies such as NGC 1569 without imposing a two-phase 
medium by hand, but still including realistic radiative cool- 
ing throughout the simulation. 

The resulting three-dimensional AMR simulations re- 
produce a number of key observational features of nearby 
starbursts, some of which have been previously captured in 
hydrodynamic simulations without gas cooling and some of 
which are unique to simulations that include supersonic tur- 
bulence. Thus, in accord with previous studies, we find that: 

• With realistic choices of star formation rates and en- 
ergy input, our simulations lead to large, bipolar outflows 
that drive substantial amounts of gas, metals, momentum, 
and energy into the intergalactic medium. These exhibit a 
"blow-out" morphology, in which the majority of the ISM 
is retained, but a large fraction of the extremely hot SN- 
driven gas escapes in a diffuse and rapidly-expanding "free 
wind." The properties of these outflows are only very weakly 
dependent on simulation resolution, and they occur without 
invoking any additional physical mechanisms such as cosmic 
ray pressure or radiation pressure on dust. 

• Most mass entrainment from the ISM occurs in the 
shear interface between the free wind and the denser ISM 
medium. Unlike other simulations with an initially homoge- 
neous medium, however, this shearing occurs both along the 
edge of the wind and from cool clouds directly in the path of 
the wind. This interacting gas leads to the majority of the 
Ha emission from the galaxy. 

• X-ray images of the galaxy show that most of the soft 
X-rays originate from the shocked material in the disk and 
from gas interactions between the hot winds and the dense 
ISM gas. Regions of X-ray emission roughly correlate with 
regions of Ha emission, but in hard X-ray images, weaker 
emission from the free wind can been seen directly, moving 
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out to large distances and escaping from the potential well 
of the galaxy. 

At the same time, in contrast with previous studies, we 
find that: 

• Structures in our simulations arise primarily from the 
interaction of shells around individual OB associations, 
which sweep up thick shells of material around more rar- 
efied pockets of hot gas. Unlike in simulations without ISM 
cooling, these dense regions persist for long times due to the 
cooling instability and the tendency for turbulence to decay 
away quickly in dense regions of gas. These effects lead to 
inhomogeneous structures throughout the starburst, which 
are far more important than the Rayleigh- Taylor instabil- 
ity in determining the outflow morphology. The result is a 
complex, chaotic Ha distribution, full of bubbles, loops and 
filaments, as observed in NGC 1569 and other outflowing 
starbursts. 

• Outflows develop not from a single large superbubble, 
but rather from the collective action of series of smaller bub- 
bles that overlap near the centre of the simulated galaxy, 
where star formation occurs most vigorously. These repeated 
outbursts open an expanding, rarefied region near the galaxy 
centre, which eats away at the denser exterior gas through 
turbulent mixing, rather than gathering it into a thin, fragile 
shell. As a result, the rarefied region drills its way outwards 
almost directly vertically, following the path along which the 
minimum amount of material separates the bubble interior 
from the intergalactic medium. 

• Blow-out occurs when the overpressured bubble regions 
from different OB associations overlap and push their way 
out into the intergalactic medium, rather than when the ma- 
terial surrounding a single superbubble becomes Rayleigh- 
Taylor unstable. This means that the strength of the outflow 
is strongly dependent on the strength of SN driving. Weaker 
outflows escape the galaxy later, are much more collimated, 
and carry far less mass, momentum, and energy, than out- 
flows from stronger starbursts. 

While each of these features are in excellent agreement 
with the Ha and X-ray morphology of starburst galaxies, 
such observations represent only a small fraction of the 
many detailed multi-wavelength constraints currently avail- 
able. In particular, a large body of emission and absorption 
line spectral data can be brought to bear in understanding 
the physics of starbursting galaxies with varying masses and 
outflow strengths. We expect future comparisons with these 
observations to not only confirm aspects of our models, but 
to lead to significant refinements. 

After all, the subgrid turbulence model presented here 
represents only a first pass at a complex and multifaceted 
problem. Our key point, then, is not that our models are 
complete, but rather that they point towards a new direc- 
tion for future research. Due to the extremely short cooling 
times in starbursting galaxies, it will be ultimately impossi- 
ble to accurately simulate them without tracking the essen- 
tial cascade of random velocities that takes place between 
the bulk motions and the thermal scale. A full understand- 
ing of galaxy outflows will only be achieved when we have 
first understood and simulated the evolution of supernova- 
driven, supersonic turbulence. 
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